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Abstract 

Supernova remnants (SNRs), as the major contributors to the galactic cosmic rays (CR), are 
believed to maintain an average CR spectrum by diffusive shock acceleration (DSA) regardless 
of the way they release CRs into the interstellar medium (ISM). However, the interaction of 
the CRs with nearby gas clouds crucially depends on the release mechanism. We call into 
question two aspects of a popular paradigm of the CR injection into the ISM, according to 
which they passively and isotropically diffuse in the prescribed magnetic fluctuations as test 
particles. First, we treat the escaping CR and the Alfvi^cen waves excited by them on an equal 
footing. Second, we adopt field aligned CR escape outside the source, where the waves become 
weak. An exact analytic self-similar solution for a CR "cloud" released by a dimmed accelerator 

strongly deviates from the test-particle result. The normalized CR partial pressure may be 

i —3 /5 

approximated as & (p,z,t) = 2 \z.\ 5 ^ 3 +z^ f 3 (p,t) exp [— z 2 /4D ISM (p) t\ , where p is the 
momentum of CR particle, and z is directed along the field. The core of the cloud expands as 
Zdif x \/Dnl (pj t and decays in time as 3? oc 2z^ (t). The diffusion coefficient Dnl is strongly 
suppressed compared to its background ISM value Dism- Dnl ~ DisM^xp (—IT) <C Dtsm for 
sufficiently high field-line-integrated CR partial pressure, IT. When IT 3> 1, the CRs drive 
Alfvi^oen waves efficiently enough to build a transport barrier « 2/\z\ -"pedestal") that 
strongly reduces the leakage. The solution has a spectral break at p = p\, r , where pt>r satisfies 
the following equation Dnl {Pbr) — z 2 /t. 



1 Also at WCI Center for Fusion Theory, NFRI, Korea 
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1. Introduction 



The generation of cosmic ra ys (CR) in superno va remnant (SNR) shocks by the diffusive shock accel- 
eration (DSA) mechanism (e.g., iDrury et a/.ll200lh is understood reasonably well up to the point of their 
escape into SNR surroundings. But making this mechanism responsible for the most of galactic CRs requires 
understanding of all stages of the CR production including their escape from the accelerators. In fact, best 
markers for "CR-proton factories" are nearby molecular clouds (MC) illuminated by prot ons leaking from 



SNRs . CRs will be visible in gamma r a ys generated by collision s with protons in the c loud lAh aronian et al. 
dl994h: lAharonian and Atovanl (fl996h: lAharonian et al\ J2008h: babici et~al\ d2009h : |Tavani et all (bold): 



Abdo et all (|2010bllah : iGiuliani et all (|201lh : Ellison and Bykovl (|201lh : iTorres et al.\ (|201 lh . Whether this 
gamma radiation is detectable with the current instruments, depends on the CR leakage rate from the source. 
The recent surge in measurements of gamma-bright SNR suggests that the sensitivity threshold have al- 
ready been crossed for at least several galactic SNRs and it is becoming increasingly timely to improve our 
understanding of the CR leakage from these objects. 

Without such improvement, it is also difficult to resolve the ongoing debates ab out the primary origin of 
gamm a-emission from some of the gamma-active remnants, e.g., RX J 1713 (e.g., Funk 2012 : Inoue et an 
20121) . In arguing for hadronic or leptonic origin, one needs to know exactly how far the CRs are spread 
from the source at a given time and with what spectrum. Indeed, strong self-confinement of accelerated CR 
may keep their flux through a remote MC below the instrument threshold, primarily (and counterintuitively) 
for powerful accelerators. Conversely, the self-confinement will enhance illumination of nearby MCs, thus 
enhancing the odds of detecting the hadronic contribution to the emission. Apart from the distance to the 
target MC, equally important is its magnetic connectivity with the CR source. Overall, the predictions for the 
emissivity of MCs near strong CR sources can differ from the test-particle results by an order of magnitude 
or more. 

To better understand the physics of CR confinement to SNRs, in this paper we consider the CR escape 
separately from their acceleration which is assumed to have faded because of the SNR age. Specifically, we 
will formulate the problem as a diffusion of a CR cloud (CRC) released from an accelerator into the ISM 
and propagating through a 'gas' of self-excited Alfven waves. At the scales larger than the initial size of the 
cloud, the solution, after an adjustment to the local environments, will become self-similar to depend only 
on the background diffusivity Dism and the integrated CRC energy (cf. Sedov-Taylor solution for the point 
explosion). 



The idea that the C RC confinement should be thought of as a self-confinement is not new (e.g JWentzel 



19741 : lAchterberglll98lh . However, the analytic solution for CR propagation uniformly valid in both the 
nearby and far zones of the initial CRC seems to be unknown. This paper is aimed at presenting and 
analyzing such solution. With some reservations, it may also be applied to a late stage of CR acceleration 
in a SNR, when the particles become disconnected from the shock and diffuse outwards. Their pressure, 
however, should still be sufficient to drive strong Alfven waves which limit their escape. We begin with a 
brief discussion of the problems with the current escape models, thus motivating the new one. 
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1.1. The need for a new CR escape model 



Traditionally, the transport of CR is treated differently inside and outside a DSA accelerator. Outside, 
CR are thought to escape as test particles with a diffusion coefficient inferable from observations. This 
picture cannot be correct near the accelerator because the CR transport must be in a self-confinement regime, 
in which the CR streaming instability diminishes their diffusion coefficient by orders of magnitude. This CR 
anisotropy instability, and particularly its nonresonant extensions, macroscopicaily driven by the CR current 
and the CR pressure gradient, have a p otential to generate magnetic field fluctuations well in excess of the 
ambie nt magnetic field, SB > B (e.g., lBeilll2004l : [Prurv and Fallelll986l : iBvkov et aZ]|2009l : llvlalkov et al. 



2010). At the very least, such strong fluctuations should justify the standard DSA assumption about the 



Bohm diffusion regime with the mean free path, (mfp) of the order of the particle gyroradius, r g , achieved 
at 8B ~ Bo. Naturally, the transport is isotropic in this regime. The ISM background turbulence, on the 
other hand, is much weaker SB 2 /Bq < 10~ 5 at the CR relevant length scales, thus resulting in the CR mfp 
> 10 5 r g (for GeV particles, or r g > 10 12 cm). It is important to emphasize that under these circumstances 
the cross-field diffusion coefficient, Kj_ is supp ressed by a l arge factor compared to the diffusion along the 



field line, K\\, i.e. K±/K\\ ~ (SB/B ) (see, e.g.. lDrurvlll983h . 



It follows then that there is a problem of describing particle transport between the self-confinement (ac- 
celerator vicinity) and the test-particle (far from accelerator) transport regimes. To circumvent this problem, 
the acceleration process has been treated separately from the particle escape usi ng one of the two devices: 



the u pper cut-off momentum and the free escape boundary (FEB) (see e.g., iReville et all 120091 : iDrury 



20111 . for recent discussions). As the names suggest, accelerated particles escape instantly upon reaching 



a prescribed boundary either in momentum or in configuration space. Their escape is a ssumed to have no 



effect on the acceleration other than through the modification of the shock struc ture (e.g., j Moskal enko et al. 
200 8). In a simplified visualization of the DSA as a 'box' process, for example (IDrury et all 19991) . the upper 



cutoff and the FEB are two sides of the box in the particle phase space. There were efforts to include self- 
generated waves into the description of particle ac celeration and escape from str ongly modified CR shocks 
( Malkov et al. 2002 : Malkov and Diamond 20061) . or in numerical treatments ( Galinsky and Shevchenkol 



2007 



201 lklFujita et a/.ll201 lh . In most approaches, however, a s udden jump in the CR diffusion coefficient 



in mo mentum space is introduced to set an upper cutoff, (e.g., 
20121) . Similar jump in coordinate space would result in a FEB. 



Ptuskin and Zirakashvilil 120051 : lYan et al. 



As long as the CR transport outside the accelerator is treated in the test particle approximation, no 
smooth transition between the CR acceleration and their escape is provided, regardless of the scenario for 
the latter. But, the diffusion coefficient rises by roughly five orders of magnitude in a transition zone that 
should be correspondingly large, unlike an infinitely thin FEB. It is this zone where the CR escape flux and 
confinement time are set by self-generated waves, thus rendering the FEB a rather implausible concept. For 
many SNRs it is then not yet possible to conclude whether the gamma-emission is leptonic or hadronic, 
should such conclusion depend on the CR escape and on the subsequent illumination of adjacent clouds. In 
a broader sense, there is a missing link in galactic CR generation between the CR acceleration (under a very 
strong self-generated wave-particle scattering) and their subsequent propagation (in a very weak interstellar 
turbulence). The goal of this paper is to establish such link. 

Our treatment below is applicable to the following two situations. In one situation, a shock accelerates 
particles continuously but some of them reach far enough or diffuse across the local field to become discon- 
nected from the shock front and have thus chances to escape. While doing so, they drive their own waves at 
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a gyroradius scale, whose amplitudes gradually decrease outwards and so does the particle density. The sec- 
ond situation is a clear-cut case for this paper that deals with a CRC released into the ISM with no ongoing 
acceleration inside the CRC. Both situations correspond to a mesoscale transport regime that is intermediate 
between the Bohm diffusion inside the accelerator and a global, quasi-isotropic CR diffusion in the Galaxy, 
where a random large scale component of the galactic magnetic field at the scale of tens of parsecs should 



wnere a random large scale component or tne galactic magnetic nela at tne scale 01 tens 01 parsecs snoula 
be dist i nguished from the regula r, spiral arm aligned toroidal field Brunetti and Codinol ( 2000l) : Strong et al. 



(|2007h : iBlasi and Amatol (|2012l ). In the mesoscale transport regime, however, both components are consid- 
ered as a large scale ambient field with a scale much larger than the CR gyroradius. Before we describe 
the next section the physical setting of the problem, it is worthwhile to emphasize some of the qualitative 
results. . 

First of all, for sufficiently high initial CR pressure (higher than magnetic pressure) and low ambient 
turbulence level, SB <C B , the spreading of the CRC strongly deviates from the usual test-particle solution. 
Instead of being controlled by only one (diffusive) scale /d ~ \/Asm^ the structure of the expanding CRC 
is more complicated. Namely, it comprises three zones, of which only the outermost has the above scaling 
but, significantly lower CR pressure. Conversely, the CR pressure in the innermost part remains strongly 
enhanced. The both zones are connected by a self-similar, 1/z CR pressure profile. 



2. CR Escape Model 



After the Sedov-Taylor (ST) expansion stage in which the shock radius increases as R s oc t 2 / 5 and 
particularly in and after the so-called pressure driven snowplow stage with a slower expansion at or below 



/? s oc t u 3 (e.g.. lBisnovatyi-Kogan and Silichl l 19951 : iTruelove and McKedl 19991) . the diffusive propagation of 
CR away from the shock becomes more important than their bulk expansion driven by the overpressured 
SNR interior. Regardless of the escap e mechanism , one of the most important parameters is the number of 
spatial dimensions involved (see, e.g.. lDruryll201ll for a recent discussion). Indeed, the 3D random walk is 
non-recurrent, so only the finite shock radius gives particle a chance to return to it. However, 2- and 1-D 
random walks are recurrent. Therefore, the first important assumption to make is, indeed, about the process 
dimensionality which immediately translates into the choice of magnetic field configuration. Note that since 
particles recede from the shock as R oc t x l 2 , they escape already at the Sedov-Taylor stage. 



Contrary to the rece nt analytic (IDrurv 1201 ll : lOhira et a/.ll201ll : lYan et a/.ll2012l) and numerical (with 
self-driven Alfven waves Fujita et al\ l201ll) treatments of the spherically symmetric particle escape, we 
consider an escape through the local flux tube. This choice is suitable for a SNR expansion in an ISM with 
a distinct large scale magnetic field direction that does not change very strongly on the SNR scale. As we 
stated in the Introduction, our goal is to fill the gap between the weak scattering propagation far away from 
the remnant and the Bohm diffusion regime near the shock, with K = K"b- The perpendicular diffusion in 
the far zone is of the order of Kj_ ~ (SB/Bq) 4 Ku <C JCn <C Kb, so that we may safely assume that in the 
unperturbed ISM the propagation is one-dimensional, along the field line. 

Closer to and inside the accelerator (or rather in the region of its past activity), K± ~ Kji ~ K"b, which 
favors the isotropy assumption. However, we may integrate (average) the equations for the CR transport and 
wave generation across the magnetic field, assuming that there is no strong expansion of the CRC in this 
directio n. The resulting problem reduces then to a l argely one dimensional propagation along the field line 
(cf. e.g. Rosner and Bodo 1996 : Ptuskin et a/.ll2008l. where a field aligned propagation from CR sources has 
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also been adopted). A cartoon example of the basic configuration is shown in FigfT] 



The second important assumption to make is about the spatial arrangement of the initial CR population. 
As shown in FigfJJ we adhere to the idea that the acceler ation is most efficient in the quasi-paralle l shock 
geometry. It can be advocated on the theo retical grounds dMalkov and Voikll 19951: [yolk et a/.ll2003h and is 
also supported by PIC simulations (e.g., Gargate and Spifkovskyll2012 ). Therefore, we may assume, for 
example, that two 'polar cusps' of accelerated particles are left behind after the acceleration has either faded 
out or entered its final stage when particles escape faster then they are replenished by the acceleration. It is 
tempting to consider the SN 1006 as a prototype of such geometry, but the similarity is physically not quite 
convi ncing, given the young age of the latter source. On the other hand, older remnants, such as W44, do 
show (|Uchiyama et a/.ll2012h a bipolar CR escape that can also be attributed to the field aligned escape. 



During earlier, more active stages of acceleration, CR presumably fill up both the downstream and 
upstream regions near the shock. Meanwhile, the contact discontinuity (CD) behind the cloud of accele rated 
CRs must have undergone the Rayleigh-Taylor instability with strong magnetic field enhancement (IGull 
19731) . The CD expansion at late evolution stages should thus act as a piston on the previously accelerated 
CRs. Note that the CR reflecting piston was already employed in numerical acceleratio n schemes, e.g., 
dBerezhkolll996h . which might, however, somewhat overestimate the maximum CR energy (|Kirk and Dendyl 
200 lh . In the post acceleration stage, and given the significant field amplification at the CD, it is reasonable 
to assume that CRs are partially coupled to the slowly expanding flow with the reflecting magnetic piston 
behind but while escaping upstream they couple to the ISM and diffuse away. 



2.1. Basic equations, initial and boundary conditions 



Perhaps the most sy stematic quasil inear derivation of equations for the coupled evolution of CR and 
Alfven waves is given by lSkillingl (1 19751) . We use them in a si mplified on e-dimensional (along the ambient 
field, z-coordinate, FigfTJ) for m equivalent to that used by, e.g., Bell ( 197 8h and, with an interesting "beyond 
quasilinear" interpretation, by iDrurvl (|1983l) 



dt 



Pcr(p) 



d k b dP C R 
dz I dz 



(1) 



±1: 

dt 



dP C R 

dz 



n. 



(2) 



Here Ca is the Alfven velocity, and the time derivative is taken along the characteristics of unstable Alfven 
waves, forward propagating with respect to the flow of speed U : 



d B , . d 

7, = T, + V +C ^T Z 



(3) 



Eq.£[]) above is essentially a well-known convection-diffusion equation, written for the dimensionless CR 
partial pressure Pcr instead of their distribution function f(p,t). We normalized it to the magnetic energy 
density pC A /2: 



-6- 



Pcr=^-^-v//, (4) 
3 pC A 

where v and p are the CR speed and momentum, and p- the plasma density. The total CR pressure is 
normalized to d\np, similarly to the wave energy density /: 



87T 



Eq.© is a wave kinetic equation in which the energy transferred to the waves is simply the diff erence 



betwe en the total work done by the particles, (U + Ca) VPcr, and the work done on the fluid, UVPcr (|Drury 



1983). This interpretation of the wave generation indicates that it operates in a maximum efficiency regime. 
A formal quasilinear derivation of this equation assumes that the particle momentum p is related to the 
wave number k by the 'sharpened ' resonance co ndition kp = eB^/c instead of the conventional cyclotron 



resonance condition kpn = eBo/c (lSkiilinpl975r) . (note that here k = k\\). We have included only the linear 
wave damping T and we will return to the possible role of nonlinear saturation effects later. We assume that 
dPcB./dz < at all times, so that only the forward propagating waves are unstable. The latter inequality is 
ensured by the formulation of initial value problem in each of the following two settings mentioned earlier 
in this section. In the first setting, the initial distribution of the CRC is symmetric with respect to z = 0, so 
we can consider their escape into the half-space z > 0. The second setting is when the cloud is limited from 
the left by a reflecting wall (CD). The appropriate boundary condition is dPc^/dz = at z = in both cases. 

Restricting our consideration to the case of coordinate-independent damping rate T, we obtain the 
following ('quasilinear') integral of the system of Equations £[]) and ©: 

Here Pcro (z) and Iq (z) are the initial distributions of the CR partial pressure and the wave energy density, 
respectively, and z' = z — (U + Ca) t. 

Using the quasi-linear integral, the system of Equations (Q]) and (O can be reduced to one nonlinear 
convection-diffusion equation for e.g., wave intensity / (z,t). We will use dimensionless variables measuring 
the distance z in units of a, which is the initial size of the CRC. The time unit is then a/C\. Note, however, 
that, as the acceleration is assumed to be inactive, the particle momentum p enters the problem only as a 
parameter but, the initial scale-height of the CRC a generally depends on p. It is also convenient to introduce 
a new variable for the wave energy density 

W = C -^I (6) 
and similar relations for Io and Wo. The equation for W then takes the following form 

dW d 1 dW d _ , , 

^o(z) (7) 



dt dzW dz dz 
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We have neglected the wave propagation along the characteristics given by eq.© and the wave damping, 
assuming that these processes are slower than the CRC diffusion: U \C^,aT <C K^/al ~ cr s /al. The function 
is defined similarly to W above 

n CR (8) 



*b(p) 

and, again, the index refers in Equation (Q to the initial CR distribution Pcro(z)- Thus, according to 
eq.©, for the dimensionless CR partial pressure we have 



This quantity is governed by the equation 



d W 

& = & - In— (9) 
oz Wo 



d&> _ d 1 d&> 

~df ~ TzWlh'' ( 
but its solution can be written down using the integral given by equation (O, after the solution to Equation 
(Q is obtained. 

While letting Ca — > in the wave-particle collective propagation, we utilize the fmiteness of Ca in 
determining the boundary condition for eq.© at z = as follows. Returning from the nonlinear diffusion 
equation eq.O to eq.©, one may see that the wave energy does not actually "diffuse" but it is generated 
locally by particles that diffuse. Therefore, even if the neglected wave diffusion had spread the waves to 
the point of their stability viz. z = 0, they would be convected away from it by virtue of Ca > 0. Recalling 
the symmetry (reflection) boundary condition dPcR/dz = at z = 0, we must set a fixed boundary value 
W (0,?) = Wb, where Wo is an initial (small) wave noise. It is instructive to investigate the case in which 
the initial noise is the same throughout the entire half-space z > (this limitation can be straightforwardly 
relaxed), so that the waves will be generated entirely by escaping particles, thus emphasizing the self- 
confinement. The second boundary condition is set by W — >■ Wo for z — > °°, and all t < oo. Note that in 
general, Wo = Wo (p). These conditions determine the boundary value problem given by eq.(|7]) completely. 
However, we are interested primarily in diffusion of CRs outside the region of their initial localization, viz. 
at z > 1, where the source term in this equation vanishes. Therefore, the problem given by eq.(|7]) can be split 
into two separate problems, one in < z < 1 and another in 1 < z < 00 domain with the following junction 
conditions 



4-lnW 

oz 



= 0. (11) 

l- 

These are the continuity conditions for both the wave energy density and pressure across the edge of the 
initial CR localization. 



i+ 

= W 

l- 



-8- 



2.2. Self-similar solution outside the region of initial CR localization 

The nonlinear boundary value problem in the region z > 1 given by Equations (Q and (fTTT ) describes 
the CR propagation outside the region of their initial localization. This problem can be solved exactly using 
the following self-similar substitution: 

W = ^w(Q, where£ = ^ (12) 

Submitting this to eq.© with = yields a = 2/3 — 1. The boundary condition at infinity (W — > Wo / 
0, z — > °°), on the other hand, requires a = 0, so that the equation for w reads 

d \ dw Cdw n 

Ji»H + 2TC° <13) 

with £ = zj\ft. It is interesting to observe that this equation has a simple special solution w = 2 /C 2 = 2t /z 2 - 



This solution describes a stationary particle distribution 3? oc \ji that is obviously related to the lBelll (11978b 
asymptotic particle distribution in self-excited waves far away from a shock. In both cases, this is a singular 
limit of the problem as it implies a zero background turbu l ence l e vel, Wp = an d thus the number of particles 
in z > half-space being infinite lLagage and Cesarskvl dl983h : iDruryl (119831) . Physically, this solution is 



different from what we will find below in that it requires a permanent source of CR at the origin so that their 
flux steadily drives waves that linearly grow in timeU As we shall see, this simple special solution is an 
important singular limit of a more general and completely regular solution. 

The general solution to eq.([T3l can be found in quadratures by swapping £ and w as dependent and 
independent variables and introducing an auxiliary function V (w) by the following substitution 

^ = -V2w 3 / 2 V(w). (14) 

"5 

The equation for V (w) takes the following form 

This equation can be easily integrated, so that the function w(Q can be subsequently found from eq.(fT4l. 
The first integral of eq.(fl3T) is as follows 

{^) 2 = R ^ sL 4 (V2 - 2 '" V -^ <16) 

where q is an integration constant. From eq.([T4l and from the boundary condition w — > Wo at £ — >■ oo we 
infer that V (Wo) = 0, so that from eq.([T6l) we find 



2 There is no convection with th e flow towa rds the origin (U = 0) in our case, which makes the solution unsteady, in contrast 
to the corresponding DSA problem ([Belli 1 9781) where the solution is steady because of the convection, but the number of particles 
upstream is still infinite in both cases. 
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fo V - 



W e V R ( V ') 

The function w (Q is then determined by the following relation for £ (V), eq.([T4l): 



w 



(17) 



c 



dV 



1 rV 

2 JO 



(18) 



Using the identity l/V = V — 2dR/dV and integrating eq.([T8l) by parts, after simple manipulations the last 
equation rewrites as follows 



W^( 2v ^ +y ) (i9) 

It is useful to introduce an auxiliary constant V\ related to the constant q and being the smaller of the two 
roots of R (V), that is R (Vi ) = 0. The requirement of a real zero of the function R (V) ensures mapping two 
'copies' of the domain of V onto the full range < £ < °°, whereby £ — > oo for V — > 0. Indeed, £ (V) in 
eq.(fl9l) diverges at V = as ~ \/—\nV. However, increasing V from V = to V = V\ does not cover the 
full range < £ < oo yet, as may be seen from eq.(fT9l: £ decreases from oo only to £i = V\ y/2/w (V\) > 0. 
To continue the integral curve to £ = 0, it is necessary to switch branches of \^R at V = V\ in eqs.(fP7T) 
and (fT9l and so continue the integral in eq.(fT71) back to V < V\ (along the second copy of the V-domain). 
Decreasing then V from V = V\ to V = Vq = exp(— q/2) brings the integral curve to the point £ = 0, 
w = w max = = 0). The explicit formal representation of the general solution of eq.([T3T> is given in 
Appendix |A] along with a numerical example. Equations (fTTT ) and ( fl9T ) determine the solution w(Q, where 
the integration constant Vq (or q) should be obtained from the matching condition, eq. (fTTT> . This can be done 
by considering eq.(|7]) inside the region of initial CR localization, i.e. < z < 1, where its r.h.s is nonzero. 



2.3. Solution inside the initial CR cloud 



Given an initial distribution 3?o{z), equation (0 can always be integrated numerically in the finite 
domain < z < 1, so that the full solution will be obtained from the boundary conditions and from the 
results of the previous section. However, according to Equation (fTOl) . already for t > Wo, where Wo <C 1, the 
initial profile (z) will be redistributed over the unity interval in such a way as to approach a quasi-steady 
state in which the flux W Q ~ l dW /dz in Equation ([7]) through z = will be balanced by the source integral, 
(0) (recall that (1) = 0). The particle flux through the z = 1 boundary decays as f -1 / 2 with time 
(since \dw/dQ < oo at £ — > 0+, see the preceding section or Equation (l23l below). Therefore, for the 
self-similar stage of the cloud relaxation we can write 




where the "integration constant" B depends slowly (in the above sense) on time. Using the first matching 
condition in Equation (fTTT) . i.e., the CR pressure continuity, we can specify B (t) as follows: B oc ? -1 / 2 — )• 
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as t — > oo (see Equation 11231 below). This determines the self-similar (outer, z > 1) solution by the second 
matching condition in Equation (fTTT ): 

w (0) = lim w (-^\ = w max = W e n (21) 
Here we have introduced the integrated partial pressure as follows 

n= f &> dz. (22) 



o 



The function B (?) and thus the internal solution W (z,t) in Equation (1201 may be obtained using this equa- 
tion, the first matching condition in Equation (fTTT ). and Equations ([T4l 



Bit) 



1 d 



lnw 



V2w/fV 



(23) 



Note that the particle pressure at z > 1 is completely determined by the turbulence level w, eq.©. 



Now we can determine the integration constant q introduced in the preceding section. From Equations 
(|2~TT ) and (1A3I) we obtain the following equation for q 



rVi rVi 

/ dV/y^R(V)+ dV/^R(yj = n, (24) 
Jo Jv 

where R (Vi ) = 0, R (Vb) = — Vb/2, while /? (V) is given by Equation (fT6l ). In the most interesting case ~ 1, 
it is convenient to use the constant e 2 = (q — 1) /2 in place of q. In this case Il3>l,e<Cl,Vi~l — £ and 

Vi l-e 

-== ^ ^2 3 / 2 ln- + ^(l), 

^(y-l) 2 - e 2 £ 

so that the turning point of the solution at V = Vi approaches the critical point V = 1 , where R has a minimum 
(Fig©: 



n = 2 



Vr~ 



For n <C 1, we can write instead, V\ ~ Vo (l + V 2 /2) where Vb = exp (—q/2). 

To conclude this section, the self-similar expansion of the CRC described by Equations (PTTfTOl) is con- 
trolled by two parameters. One parameter is the background turbulence level Wo (p) in the media into which 
the cloud expands. The second parameter is the integrated pressure of the cloud, IT. Although the initial 
wave energy density inside the cloud (z < 1) is likely to be higher than Wq, we have adopted the background 
value Wq also inside the cloud for simplicity, as it should not influence the self-similar CR propagation out- 
side the cloud. In addition, efficent wave generation by a dense expanding CRC, renders the initial value for 
the wave energy density inside the CRC unimportant. 
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3. Analysis of the solution 

Once we have the solutions outside and inside the region of initial CR localization, we may precisely 
calculate how fast particles escape from this region. Two convenient characteristics of the escape process are 
the half-life time and the width of the CR distribution. We begin our analysis of the solution from computing 
these simple quantities and turn to the details of the CR escape afterwards. 



3.1. Half-life time and the width of the CR cloud 

The most concise characterization of the escape may be given by looking at the following two parts of 
the (conserved) integral particle pressure: 

n = n + rii = const 

where 



1 oo 

n = j @>dz and Ui = j 



8*>dz 



refer to the regions inside and outside of the initial CRC, respectively. Recall that at t = 0, Ilo = IT and 
111 = 0, while at t = oo, an opposite CR distribution is reached: Ilo = and IIi = IT. Note that Equation 
(|2TT) specifies the total work ultimately done by particles on waves, while they diffuse from < z < 1 to 
1 < z < oo. As n is conserved, it is natural to define the CR confinement time as the time at which Ilo (Oi ) 
drops (raises) to a half of IT. Substituting W from Equation (fT2l) into Equation ©, we obtain 

For the half-life time t = t\/2 we may use the following equations: 

ri! (f 1/2 ) = n ( , (* 1/2 ) = ^n. (25) 

In the simple test particle case (n <^ 1), the half-life amounts to (see Appendix IBT) 

1 

4c 

with a w 0.48. In the opposite case n 3> 1, for ti/ 2 we obtain 



h/2 ~ -^Wq (26) 



, 2 _ W\ W\ 

(see Equation |fT9l ). Note that for IT 3> 1, Vi pa V\ii, (see Appendix |B]) and for the half-life time we obtain 
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h/2 = \w e u / 2 (28) 

We may extend the last formula with a 5 percent accuracy also to the linear case (IT <C 1), given by Equation 
(l26l ). It is clear that the nonlinear delay factor exp (IT/2) slows down the escape considerably, compared to 
the test particle solution. 

Another important characteristics of the particle escape is the spatial width of their distribution. For- 
mally, the self-similar solution is scale-invariant ~ 2/£, for n S> 1 over the most part of the spatial distri- 
bution of the partial pressure (see below). Therefore, to characterize the width we use the point £in, related 
to the half time fy 2 as follows: 




In physical coordinate z, this point moves outwards as z.\/2 = Ci/2V^ starting from z = 1 at t = t\j2- The 
expansion rate of a dense CRC (IT » 1) is thus exponentially low. 



3.2. Spatial distribution of the CR cloud 

Considering the spatial distribution of the spreading CRC in detail, it is convenient to start with the 
region far away from the source, where W — > Wq. This asymptotic behavior corresponds to the case V <C 1. 
Introducing an auxiliary function U (Q 



U 



-t— rlnw = V2wV 



(29) 



that is related to the particle pressure through & = U (Q / \ft, and using Equations (fT6l ) and (fT9l , we obtain 
for U the following expression 



U (Q = y/2WoV e 



WoC 2 /4 



(30) 



Therefore, the asymptotic CR pressure depends on the following two parameters: the ISM background 
diffusivity Wq (p) = W~ l (p,z = °°) (see Equation ifTOl ) and on the CR source pressure II (p), but only 
through the parameter Vq = exp (—q/2). The latter grows linearly with IT <C 1 but it saturates when II 3> 1: 



Vb 



exp 



7b n > 

I _4 e -n/V5 



n«i 
n> i 



(31) 



For practical calculations, it is more convenient to combine both asymptotic regimes into the following 
simple interpolation formula: 



Vo 



^/nf 2 + (y-e + \e-^ 



5/2 



-2/5 



(32) 
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which not only recovers both IT <C 1 and IT 3> 1 regimes but also reproduces the transition zone accurately. 
The latter property is achieved by the choice of numerical parameters of the interpolation, i.e., 5/2, 3/2 and 
2/5. The quality of the interpolation may be seen from FigJH where the formula in Equation (l32l is shown 
against exact numerical points calculated from Equation (l24l) . The partial pressure (z,t) is then given by 



& = u/y/i = y/mjtv (n) e - w ^l At 



(33) 



Note, that for strong sources (11 3> 1) the CR density at large distances becomes independent of the source 
strength. This is the regime of an efficient ablative escape reduction in which a small (and IT-independent!) 
number of leaking particles leave behind enough Alfven fluctuations to limit the leakage of particles remain- 
ing in the source to exactly the rate given by Equation (|33T >. 

The spatial distribution of the CRs becomes even more universal closer to the origin where it falls off 
as & « 2/z before it turns into an innermost flat-top part of the entire distribution for £ 2 < Dnl (/?)• We 
have introduced the "self-confinement" diffusion coefficient Dnl as 



D 



NL 



777 £>iSM<? 
V 



-n 



V 2 w n 



(34) 



with Dism = Wq i . To obtain this intermediate (Dnl < C < Asm) part of the CR spatial distribution we 
expand the analytic solution given by Equations (flTT ) and (fT9l > in small V\ — V <C 1 . Using the expression 
for the self-similar CR pressure U from Equation d29l ) we obtain for this quantity the following expansion 
nearC = Ci = C(Vi): 



2 

U = X 



x/2 



1/V2 



%V2 



+ .. 



(35) 



where w\ = 2Vf /Q. Note that while the solution has a branching point in auxiliary variable V at V = Vi, 
it is regular and single- valued function of the physical variable £, throughout the entire half space £ > 0, 
including the branching point of C(V) at V = V\, as it should be. We also notice that while £ decreases 
starting from £i, the solution grows slower than l/£, leveling off towards the origin. The above expansion, 
however, becomes inaccurate for smaller £ and we alter it below. 

Now we turn to this innermost part of the distribution where the particle diffusion coefficient is most 
strongly decreased. It is convenient to expand the solution in a series in t, = V /Vo -1 < 1. Using the 
representation of w given by Equation (IA2b . we obtain 



where a ■ 
Equation 



w ■ 



b — a 



-2b 



1 - V 2 ) / (l + Vq) and b = Vo/J 1 + V 2 . Then, using the expression (lA4"b for small £ and 
, we obtain for U (Q = &\ft the following simple result 



U 



V2w max V 
l + w max £ 2 /4 ! 



(36) 
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valid for w max £ 2 ^ 1- Note that, with the same accuracy, the denominator can be recast in the standard "diffu- 
sion" form, exp (w max £ 2 /4), in which the CR diffusivity is diminished by a factor Wb/w max = exp (—IT) <C 1 
compared to its background level. If, however, w max > oCa/ Kb (see Equation J6l), a better approximation 
would be to simply replace the diffusion coefficient by its Bohm value, as the neglected nonlinear wave 
interactions (see Secj4j) should render 8B ~ Bq. The result in Equation (l36l ) should then be replaced by 



The overall distribution of U (Q = &^Jt is well reproduced by the following parametrization, Figj5] 



where Vq (n)is given by the interpolation in Equation (l32l ). The last formula recovers the limiting cases 
given by Equations (l30l ). (1351) and (1361 ) except a slight modification of Equation (1301 ) at large £ , so that an 
extra 1 / ^-factor is accepted in the interests of the simplicity of the parametrization. According to FigJU 
however, this deviation from the correct asymptotic expansion of the exponential tail of the distribution is 
not really noticeable. 

The overall profile of the partial pressure (Fig© unveils the CR escape as a highly structured process. 
It comprises a nearly flat-top core, a 1/^ pedestal and an exponential foot. We already noted that the escape 
rate of CR in the foot saturates with the source strength for II S> 1 . This is easy to understand as the foot 
is separated from the core (source) by the pedestal, where the CR transport is self-regulated so that a fixed 
CR flux streams through it to the foot. The solution is shown in FigfS] as a function of £ using the self- 
similar variables U = &\ft and £ = z/\/i. Remarkably, the pedestal portion of the profile does not change 
with time also in the physical variables, & » 2/z. The core, however, sinks in as & oc l/y/i but in 

addition, it expands in z as a/Z5nI7> to conserve the integrated pressure contained in the core at the expense 
of the pedestal that shrinks accordingly. The latter disappears completely at t ~ Dism/Aml, after which the 
further escape proceeds in a test particle regime but with a smaller diffusion coefficient in the core-pedestal 
region, since waves persist even after the most of particles have escaped. 

As the test particle approximation is widely used in calculating the CR escape from their sources, we 
also show, for comparison, one such example in Fig|5] We plot the following simple test particle solution 



that can be obtained from Equation (l38l) or (l33l) with Vb (IT) taken from Equation (f3TT > or (l32l for IT <C 1. 
However, we formally extend this linear result to large IT, as this is normally done within test particle 
approaches. Such an extension clearly underestimates the nonlinear level of the CR pressure by a factor 
~ IT 1 exp (11/2) S> 1 in the core and in the most of the pedestal. By virtue of lacking self-regulation, it 
overestimates the pressure in the exponential foot by a factor IT » 1 . 



U = ^/2aC A /K B V (IT) e - aC ^ 2 ^' 



(37) 



U = 2 C 5/3 + (D NL f 6 



3/5 



e -WoC 2 /4 



(38) 




3.3. Control parameters and predicted CR flux in physical units 



The integrated partial pressure IT is the most important parameter that regulates the CR escape. There- 
fore, we consider it in some detail, returning to the dimensionful variables and rewriting IT as follows (see 
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Equations O and EH ): 



n^^^M (39) 



where is the initial size of the CRC, r g is the CR gyroradius and Pcr(p) is their average partial 
pressure inside the cloud. Although the ratio of Alfven speed to speed of light, Ca/c is typically very small 
(~ 10~ 5 — 10~ 4 ), the remaining two ratios in the last equation may be fairly large. Indeed, a/r g should 
amount to several c/Uest 3> 1> with Uest being the shock speed at the end of the Sedov-Taylor stage, as 
a scales at least as the precursor size. In fact, a significant part of the shock downstream region filled with 
CR, that have been convected there over the entire history of CR production, may significantly contribute 
to, if not dominate, the total length of the cloud a. This is p articular l y rele vant to the CR release during 



the decreasing maximum energy (reverse acceleration phase) iGabicil (120111) . In the forward acceleration 
regime (growing maximum energy) Bohm diffusion, the downstream CR scale height is similar to that of 
the upstream. A reasonable estimate for Uest is c/Uest ^ 10 3 . Finally, the CR pressure in the source should 
considerably exceed the magnetic pressure as both quantities are roughly in equipartition in the background 
ISM. Alternatively, as is usually assumed, the accelerated electrons are in equipartition with the magnetic 
energy inside the source, since they should have lost their excessive energy via synchrotron radiation. As 
electrons are thought to be involved in the acceleration at ~ 10 of the proton level, the last ratio in Equation 
(|39l is at least > 10 2 , thus raising the pressure parameter to IT ^> 1. 



One of a few other ways the parameter IT can be looked at is to express it through the average ac- 
celeration efficiency $ = 2Pcr/ pU 2 , where U is the shock velocity, appropriately averaged over the CR 
acceleration history (one may expect U > Uest)- Then IT ~ 3 (U 2 /cC A ) {a/r g ) S ~ 3A (0 2 /U E stCa) <§ , 
where we have introduced a factor A = (a/r g ) Uest/c > 1. In this form, the estimate indeed boils down to 
the acceleration efficiency S 1 with the remaining quantities (A,U and Uest) being more accessible to spe- 
cific models. Since the acceleration efficiency $ is believed to be at least > 0. 1 for productive SNRs, we 
conclude that the control parameter IT is rather large than small. The caveat here is that it might, in some 
cases, be too large to limit the applicability of the above treatment. We return to this issue in the Discussion 
section. 

Once the major control parameter is known, we can calculate the distribution function /cr of the CR, 
released into the ISM, depending on the distance from the source z and time, using the parametrization in 
Equation (l38l) . It is convenient to rewrite it in the form of the CR partial pressure Pgr in physical units as 
follows 

fe^ 08 /«Aiov . _L_f(H^Y ,2 u{Q m 

B 2 /Sn \A t lcm- 3 lGeVy V B J 

Here Ma = Uest /C a, n is the plasma number density, E is the particle energy and B is the magnetic field. 
The self-similar coordinate £ can be represented, in turn, in the following way 

/^1^^10VsV /2 A ( 4i) 
b V A E ljUG t J pc 

Once this quantity is calculated for given z and t, the CR partial pressure in Equation (1401) can be obtained 
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from Equation [38] or from FigfSJ Note that the scale of the CR pressure is determined by the maximum of 
U = U (0) = ^/2WbVbexp (IT/2), where Vq is given by Equation d32l ). The CR partial pressure also depends 
on the energy through Wo and n, apart from the factor yE in Equation (l40l) . which will be discussed 
below. Finally, the background diffusion coefficient of CR in physical units is related to the dimensionless 



diffusivity W 1 as follows 



77 1(T 4 a B ( lcm- 3 \ 1/2 cm 2 

D ISM = 6.6- 10 27 -— — -— 

Wo lpc ljitG \ n J 



sec 



3.4. The Form of the Spectrum of escaping CR 

As the transition from the flat-top part of the normalized (Equation JH) partial pressure & (p) to 
the pedestal is described by a function 3? m 2|z 5 / 3 + [Dnl 5 ^ 6 j , the pressure 3?{p) is almost 
/^-independent (the coiTesponding particle momentum distribution / CR oc K B p~ 4 /z, in the physical units, 
z = az) for such momenta where Dnl {p) < z 2 /t, for the fixed dimensionless observation point z and time t. 
At p = pt,r> where Dnl {Pbr) = z 2 /t, the spectrum incurs a break. If we approximate Dnl x P S near p ~ pt>r, 
the break will have an index 8/2, as may be seen from the above formula for The index 8 thus derives 
from the momentum dependence of the Dism (p) and from that of the coordinate integrated CR pressure 
n(p) (through the factor exp(— IT) in equation 11341 ). Furthermore, if we represent exp(— n) oc p~ a and 
Dism x P X at p ~ p^ r , so that 8 = X — a, then 3? is flat at p < pb r for 8 > and steepens to p~ 5 l 2 atp = pb r . 
Conversely, if 8 < 0, & raises with p as p~ s l 2 at p < p\, r and it levels off at p > p\, r . Note, however, that 
n is momentum independent in the important case of the p~ 4 distribution of the initial CRC with the scale 
height a(p) K"b- This case corresponds to the test-particle DSA spectrum, oc p~ 4 with a scale height a 
estimated from the shock precursor size. The break has then an index A/2 and it is entirely due to the Dism 
momentum dependence. 



4. Comparison with other approaches 



Given a number of different approaches to the CR escape, we extend our brief discussion in the Intro- 
duction to compare them with the above-developed treatment. Most of the approaches can be categorized 
into the following three kinds. First of all, a sim ple test particle approach (TP ) is feasible for a rarefied 



CRC , when wave generation is negligible (e.g., lAharonian and Atoyanl 1 19961 : iDruryl 1201 ll : lOhira et al. 



20111 ). It solves the linear diffusion equation for the CRs with a diffusion coefficient determined by a 
given (e.g., background) turbulence. Next, there are extensions that also estimate the wave generation by es- 
caping particles and determine the nonlinearly satu r ated waves, mos tly by considering MHD cascades (e.g., 
Ptuskin and Zirakashvilil 120051 : IPtuskin et a/.ll2008l : lYan et a/.ll2012h . The wave intensity, thus obtained, is 
then used to calculate the particle diffusion coefficient and their distribution. We tentatively call these ap- 
proaches modified test particle, as they still do not treat particles and waves on an equal footing. By the 
conventional plasma physics methodology this improvement of the TP models overjumps the quasi-linear 
(QL) treatment, where wave stabilization occurs through their backreaction on the particles (suppression 
of instability by pitch-angle isotropization) and not through the nonlinear saturation. As this is the lowest 
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order process in the wave energy that limits the wave growth, it is imperative to consider it first. We thus 
pursue a conventiona l plasma physics approach, s imilar to a classical problem of expansion of a hot electron 



cloud into a plasma (IRyutov and Sagdeevlll970l : llvanov et al\\\91w . As we also noted earlier, there is an 



extended region (i.e. "pedestal" self-confinement region) where the level of turbulence strongly exceeds 
its background value but is stilll not high enough for the nonlinear effects to dominate over the quasilinear 
ones. As our results show, this is dynamically the most important domain that determines the spectrum of 
escaping CRs. Simply put, the dynamics is dominated by wave generation and self-regulated particle escape, 
rendering the nonlinear wave dynamics and MHD cascades less important. It should be noted, however, that 
the instability driving pitch-angle anisotropy is supported by the spatial inhomogeneity of the CRC, so that 
the full stabilization occurs (under negligible damping) only when particles spread to infinity. 

The above-described approaches to the CR escape are charted in Figj6] The key ingredients are the 
CR particles and Alfven waves while the relevant physical phenomena are the wave-particle and wave-wave 
interactions. Under the latter we loosely understand both weakly-turbulent parametric processes (such as 
wave decay) and turbulent cascades. The wave-particle interactions comprise the unstable growth of Alfven 
waves together with their back-reaction on the CRs. Wave-wave interactions need to be included in the 
quasi-linear treatment, particularly in the core of the CRC when w max > 1 , but the exact analytic solution 
obtained in this paper will hardly be possible to extend to the fully nonlinear treatment. Therefore, at present 
we may employ the Bohm diffusion coefficient in this region, i.e. where z 2 /t < Dnl by setting Dnl ~ Kb 
and thus obtain the expression in Equation (I37T ). which formally amounts to the test-particle result with the 
diffusion coefficient decreased to its minimum, i.e. Bohm value. 

To summarize our results, we have considered the self-consistent relaxation of a CR cloud (CRC) 
injected into a magnetized plasma (ISM) under the assumption of an initially weak background turbulence, 
(SB/Bq) 2 1, so that the cross-field diffusion is negligible, Kj_ <C fCii outside the cloud and the particles 
escape largely along the field, i.e., in z-direction. The principal parameter that regulates the CR escape from 
the cloud is identified to be their coordinate-integrated partial pressure II, given e.g., in Equation ( |39l ). 
Resonant waves of the length ~ r g 1 (p), driven by the run-away cloud particles, are found to confine the 
core CRs very efficiently when IT 3> 1 - 

The resulting normalized (Equation JH) CR partial pressure profile @* comprises the following three 
zones: (i) a quasi-plateau (core) at small z/y/i < VAml of a height ~ (Dnl*) 2 > which is elevated by 
a factor ~ n _1 exp(II/2) 3> 1, compared to the test particle solution because of the strong quasi-linear 
suppression of the CR diffusion coefficient with respect to its background (test particle) value Disivl Dnl ~ 
DisM^xp (—IT) (ii) next to the core, where \/D^l <z/\/t < \/Djsm, the profile is scale invariant, & « 2/z. 
The CR distribution in this "pedestal" region is fully self-regulated, independent of IT and Dism for IT ^ 1 , 
(iii) the tail of the distribution at zj\[t > \/Djsm is similar in shape to the test particle solution but it saturates 
with IT 3> 1, so that the CR partial pressure is <=c (DjsmJ) 2 exp ( — z 2 /4Dism?), independent of the strength 
of the CR source IT, in contrast to the test-particle result that scales as oc n. Because of the CR diffusivity 
reduction, the CRC half-life is increased and the cloud width is decreased, compared to the test particle 
solution. 



Depending on the functions IT (p) and Dism (p), the resulting CR spectrum generally develops a spec- 
tral break for the fixed values of z and t such that z 2 jt ~ Dnl {p) ~ DisM^xp (—IT). 



-18- 



5. Discussion and Outlook 



The CR escape from both active and fading accel erators (old SNR) is being actively studied t hrough di- 
rect observations of CR illuminated molecular clo uds (|Aharonian et al . 2008 ; Abd o et al Giuliani et al. 
2011uAleksic et a/.ll2012l : IUchiyama et a/.ll2012r) . To date, most of the information is obtained from the old 
remnants and they consistently show a broad spectrum CR escape. This is clearly at odds with an intuitive 
high energy biased escape, seemingly justified by the higher mobility of energetic CRs. Indeed, as CR diffu- 
sion coefficient grows with momentum, the test-particle solution predicts a low-energy cutoff to be present 
due to the factor exp [— z 2 /4Dism {p)t\ in the CR distribution at a certain distance from the source. In a 
combination with a steep power-law or a favorably placed upper cutoff, the escape flux narrowly accumu- 
lates towards the maxim um energy. The momentum dependent CR mobility underlies most of the curren t 



CR escape models (e.g., iPtuskin and Zirakashvilill2005l : 



Zirakashvili and Ptuskinll2008l : lGabici et aZ.ll2009h . 



Although the same exponential factor is present in the QL solution obtained in this paper, it pertains 
only to the farthermost zone, where the CR partial pressure is much lower (by factor II 1 <C 1) than the 
test particle prediction for the same distance from the source and time elapsed from the CR release. Closer 
to the accelerator, in an extended scale-invariant zone where the CR level is much higher and decays as 
slowly as l/z, the escape mechanism is different from the one controlled by the mere energy dependence 
of the CR diffusivity. It is self-regulated in such a way that, if particles leak excessively in some energy 
range, they also drive stronger resonant waves to reduce their leakage and vice versa. As a result, the 
overall escape spectrum relaxes roughly to an equipartition of the CR partial pressure in momentum, e.g., 
fcR x P (with important deviations described in Sec l3.4l) . which also balances the driver (gradient of 
CR partial pressure) with the generated waves. No low-energy leakage suppression therefore occurs. The 
fundamental difference of this leakage mechanism from the test particle one is that it is entirely controlled by 
self-generated rather than prescribed waves, or by waves derived from other energy sources, such as ambient 
MHD cascades. That is why it predicts energetically much bro ader leakage than many other approaches do 
(e.g., IPtuskin and Zirakashvilill2005l : Ellison and BykowOllI) 



It should be admitted that the self-confinement solution obtained in this paper is strictly valid for a 
stopped accelerator, so that there is no CR energy growth and associated strong plasma flows, such as those 
found near shocks. Therefore, care should be exercised in comparing this solution with the standard DSA 
predictions. Nevertheless, particles that escape from the DSA through free escape boundary (FEB) or the 
upper momentum cutoff, are likely to propagate further out diffusively. Their escape is, in fact, enforced by 
imposing an ad hoc sudden jump in the diffusion coefficient at a specific momentum (or FEB position) that 
advances in time according to the acceleration rate and shock evolution. Despite this jump, the CR phase 
space density just at the cutoff (or FEB) must be continuous^ i.e., still high enough to drive waves while the 
CR escape. From this point on, the solution obtained in this paper may be applied and compared with the 
DSA predictions, as the waves are driven locally both in momentum and in coordinate space, and the only 
relevant requirement is that particles do not interact with the shock, as implied in most escape models in the 
first place. But, the feedback from the self-generated waves on the CR escape is not included in the models 



3 This state ment is strictly valid for a FEB imposed by enhanced diffusion in coordinate space. In the case of a jump in mo- 
mentum (e.g.. IPtuskin and Zirakashvilill2005h . the situation is more complicated in that the particle distribution should become 
increasingly anisotropic in pitch angle, as the escape is assumed one-sided (upstream). However, an introduction of weak Fermi-II 
acceleration (diffusion in momentum) would validate this statement also in this case. 
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that predict a peaked energy escape. 



Furthermore, the self-regulated escape solution shows a gradual increase of the CR diffusion from a 
low (Bohm) to high (TP) regime, and most of this growth occurs where the solution is scale-invariant, i.e. 
the partial pressure 8? oc \/ z and the CR diffusion coefficient D NL °= z 2 , thus keeping the flux across the scale 
invariant (pedestal) region —D^iV^ 1 ps const. Both scalings are clearly inconsistent with a sudden jump in 
D NL with the corresponding jump in Moreover, collapsing the scale-invariant region to a point would 
not only change the CR escape flux considerably but, an extended and observationally important region of 
the enhanced CR pressure would be completely lost (see Fig[5j). This region (which we loosely dubbed 
"pedestal" by analogy with the improved confinement regimes in magnetic fusion devices, primarily in 
tokamaks e.g., Wagner et al. 1982 : Hinton and Staeblerlll993 : Diamond et al. 1995 : Malkov and Diamond 
20081 : iMaggfeoid) may be detectable when it overlaps with molecular clouds. According to the test-particle 
theory the CR density decays as t~ d l 2 in the region of their initial release, with d being the dimensionality 
of the escape (d = 1 for escape along the field). In the self-regulated escape, the CR density stays constant 
in time in the pedestal region, where Pgr °= l/z, until this region is overwhelmed by the expanding central 
plateau where Pqr decays as 1 / y/t. Besides that the CR density is significantly higher there than in the 
test-particle case, FigfS] which is the prediction that may soon become testable. 

Recent detailed observations of the SNR W44, W51C, IC 443 and W28 surrounded by dense MC 
provide good examples to study possible CR escape scenarios ( Aharonian et a/.ll2008l : Abdo et 0020 10b 3fi 
Giuliani et a/.ll201 ll : lAleksic et a/.ll2012l : lUchiyama et a/.ll2012l) . First of all they invariably show spectral 
breaks that, howe ver, may be understood in terms of the interaction of accelerated protons with a partially 
ionized dense gas (|Malkov et al 2005, 1201 ll) . The indices below and above the breaks are consistent with the 
following two scenarios. Namely, CR protons may reach the molecular cloud while still being accelerated 
at a SNR shock, or they may escape with a similar spectrum, p~ 4 . As we have shown in the present 
paper, such escape occurs via the CR interaction with self-generated waves. Another important aspect of 
these observations regards the morphology of the interaction. Of particular interest is the recent analysis 
of Fermi-LAT results r evealing two bright spots of gamma emission adjacent to the central source in W28 
(|Uchiyama et a/.ll2012r) . Their distinct bi-polar appearance may be indicative of a CR escape along the local 
magnetic field. 

Note that the anisotropic diffusion of cosmic rays in the form of bipolar CRC may result in quite specific 
morphologies of extended gamma-ray images - the imprints of cosmic rays interacting with the surrounding 
diffuse gas are generally concentrated in dense molecular clouds. Since the gamma-ray flux is proportional 
to the product of densities of CRs and the diffuse gas, we should expect a rather general correlation between 
the gamma-ray fluxes and the column densities of the interstellar gas. However, it is obvious that in the 
case of propagation of bipolar CR clouds through an inhomogeneous clumpy gaseous environment, the 
gamma-ray intensity contours can significantly deviate from the CO and 21cm maps characterizing the 
spatial distributions of the molecular and atomic gases, respectively. At the same time, the brightest parts of 
the spatial distribution of gamma-rays should correspond to the regions where the CR cloud overlaps with 
dense gas clouds inside the magnetic flux tube connected with the CR source. 
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A. Details of self-similar solution 

The full expressions for the solution w(Q represented in a short form by eqs.([T71) and (fT8|) can be 
written as follows 



w = Wo\ V] 



e 5Xdv/VW) 



C>Ci 

In some cases, it is convenient to represent w in the domain < £ < £1 as 



(Al) 



w = w max e Jv o 
where Vq = exp (—q/2), i.e., £ (Vb) = 0, and 



-Sldv/y/W) 



(A2) 



W m ax = Woe 



lo'dV/y/m+f^dV/y/W) 



For the coordinate £ we have 



2/ C > Ci 

v-2^0 7 ) o < C < Ci 



(A3) 



(A4) 



Throughout this Appendix, the positive branch of VR is used. The partial CR pressure can be represented 
according to eq.® as follows 



where V (Q is defined by eq. (|A4|) . 



B. Half-life of the CR cloud 

Equation (l25l) can be continued as 

n 1 (^ /2 ) = in = iln^i^ln^ (Bl) 

where we have introduced a "half-life" amplitude Win, which can also be associated with a point V1/2 (see 
EquationlEl) 



and with the corresponding self-similar coordinate = 1 / y/h/2- m other words, W1/2 = w (C1/2) = 
w (V1/2). Note that W\/2 can be represented as a geometric mean W1/2 = vWow max (cf. equation |[24l ). so 
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that in terms of the function lnw(Q, the point £]/ 2 corresponds to the FWHM. The quantity V\/ 2 can thus 
be obtained from the following equation 

2 / dV/y^R(V)= / dV/VW)+ / dV/^/RAV) (B2) 
Jo Jo Jv 

In the case of a weak CRC, II 1, we obviously have Vin < Vq ~ V\ <C 1, so that substituting /? from 
equation (fTBT ). we have 

erfcflnJ--^=i 
V Vi /2 2; 2 

or 

where erfc(a) = 1/2, so that a « 0.48. Substituting then V1/2 into £1/2 from equation (IA4b we obtain the 
result given by equation (l26l) . 

In the case IT » 1, the integrals on the r.h.s. of Equation (IB2I ) are dominated by the upper limit, so the 
integral on the l.h.s. must also be and we deduce V\/ 2 ~ Vi, from which we obtain the nonlinear CR half-life 
in Equation (1281) . 
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Fig. 1.— 



CR escape along the magnetic field Bq from the two polar cusps of SNR with a stalled blast wave. 
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Fig. 2. — Function R(V) and the limits of integration in equations (fTTT ) and (|A lb depicted for e = 0.3. 
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Fig. 3. — Analytic vs numerical solution of eq. (fl4l) . The gap in the analytical curve encloses the branching 
point of the solution at V = V\, eq. dATI ). w = 0.19, w m = 0.9 
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Fig. 4. — Squares: Function Vb (II), obtained from Equation (|24l . Line: interpolation given by Equation 
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Fig. 5. — Spatial distribution of CR partial pressure (as a function of £ = z/y/t, multiplied by y/i) shown for 
integrated values of this quantity IT = 3.6; 6.7; 10.1 and for for the background wave amplitude Wo = 10~ 4 . 
Exact analytic solutions are shown with the solid lines while the interpolations given by Equation (I38T ) are 
shown with the dashed lines. For comparison, a formal linear solution for n = 10.1 is also shown with the 
dot-dashed line. Note the three characteristic zones of the CR confinement: the innermost flat top core, the 
scale invariant (I/O pedestal, and the exponential decay zone. 
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Fig. 6. — Three different approximations for studying injection of accelerated CRs into ISM: Test particle 
(TP), modified TP (with unstable wave growth and nonlinear wave evolution), and quasi-linear (QL) (with 
self-consistent time dependent wave-particle interactions). 
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